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Abstract. In this colloquia review we discuss methods for thermal transport calculations for nanojunctions 
connected to two semi-infinite leads served as heat-baths. Our emphases are on fundamental quantum 
theory and atomistic models. We begin with an introduction of the Landauer formula for ballistic thermal 
transport and give its derivation from scattering wave point of view. Several methods (scattering boundary 
condition, mode-matching, Piccard and Caroli formulas) of calculating the phonon transmission coefficients 
are given. The nonequilibrium Green's function (NEGF) method is reviewed and the Caroli formula is 
derived. We also give iterative methods and an algorithm based on a generalized eigenvalue problem 
for the calculation of surface Green's functions, which are starting point for an NEGF calculation. A 
systematic exposition for the NEGF method is presented, starting with the fundamental definitions of 
the Green's functions, and ending with equations of motion for the contour ordered Green's functions 
and Feynman diagrammatic expansion. In the later part, we discuss the treatments of nonlinear effects in 
heat conduction, including a phenomenological expression for the transmission, NEGF for phonon-phonon 
interactions, molecular dynamics (generalized Langevin) with quantum heat-baths, and electron-phonon 
interactions. Some new results are also shown. We also briefly review the experimental status of the thermal 
transport measurements in nanostructures. 

PACS. 05.60.Gg quantum transport - 44.10.+i heat conduction - 65.80.+n thermal properties of small 
particles, nanocrystals, and nanotubes 



1 Introduction 

The first quantitative description of the phenomenon of 
heat conduction was given by Fourier in the early 1800s, 
which states that the heat current is proportional to the 
temperature gradient, J = -kVT. The coefficient, k, is 
known as the thermal conductivity. Together with energy 
conservation, the temperature field of a macroscopic body 
satisfies a diffusion equation. Debye proposed simple ki- 
netic theory to express the thermal conductivity in terms 
of a product of specific heat, velocity, and mean free path 
of the phonons. A more complete theory is given by Peierls 
PP, based on a Boltzmann equation for the phonons. The 
Boltzmann-Peierls equation approach for understanding 
the thermal transport has been a standard approach, and 
has been perfected by many people [2]. There is also re- 
cent work for a more rigorous derivation of the Boltzmann 
equations [3]- 

In recent years, much research has been done on meso- 
scopic to microscopic systems, including thermal trans- 
port |4|5|6j . for the obvious reason of relevance to the 
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miniaturization in the electronic industry. At small scale, 
some of the concepts familiar at macroscopic regime may 
not be applicable, e.g., the concept of a distribution func- 
tion of both coordinate and momentum used in the Boltz- 
mann equation. At the nanoscale, where atomic detail 
becomes important, we expect to replace the Boltzmann 
equation by something more fundamental. One possibil- 
ity is quantum Boltzmann equations [7|8j . Another class 
of methods for thermal transport problems starting with 
atomic models is molecular dynamics (MD) 9 10 11 . The 
thermal conductivity can be computed in equilibrium based 
on the Green-Kubo formula, or through a direct computa- 
tion of the thermal current in a nonequilibrium setting. A 
variety of systems has been simulated with MD for nanos- 
tructures and interfaces. However, since MD is purely a 
classical method, quantum effect can not be taken into 
account. 

In this review, we consider a bottom-up approach where 
we start with a lattice model with only the harmonic inter- 
actions. Such models explain well the thermal transport 
properties when systems under consideration are smaller 
than their mean free path of the phonons. This regime 
is known as ballistic heat transport regime. This will be 
the contents of Sec. [2] In this section, we discuss the con- 
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cept of universal thermal conductance for perfect ballis- 
tic systems, and address the issue of computing the ther- 
mal current if the system is not a perfect periodic lattice. 
The result is given by the Landauer formula. Effective 
computational methods will be reviewed for the trans- 
mission coefficient used in the Landauer formula. What 
will happen if we introduce interactions for the phonons? 
We discuss this issue in Sec. [3] The problem of nonlinear 
interaction can be solved, at least in principle, with the 
nonequilibrium Green's function (NEGF) method, which 
is more fundamental than the Peierls approach, and with- 
out conceptual difficulty. In practice, the NEGF method is 
computationally extremely intensive if the nonlinear prob- 
lem is treated by mean field approximations. We also re- 
view a new molecular dynamics method where quantum 
effect can be taken into account approximately. In par- 
ticular, the ballistic regime can be simulated correctly. 
Electron-phonon interaction, as well as recent experimen- 
tal progress, will be reviewed briefly. 



1.1 Models in this review 

For notational consistency, in this review, we'll discuss a 
model of junction connected to two leads using the follow- 
ing Hamiltonian [12ll3j . The model is in fact very general. 
In order to have a form for the Green's function similar to 
that of electrons, we use a transformation for the coordi- 



nates, Uj 



where x n 



is the relative displacement 



of j-th degree of freedom having the dimension of length, 
and we call Uj the mass-normalized displacement. In this 
way, the kinetic energy is always of the form ^u T u (where 
the superscript t stands for matrix transpose), and the 
mass information is transferred completely to the spring 
constants. We introduce a superscript a to indicate the re- 
gion. Then is the deviation from equilibrium position 
for the j-th degree of freedom in the region a; a = L, C, i?, 
for the left, center, and right regions, respectively. The 
quantum or classical Hamiltonian is given by 



C\TTrCR n ,R 



a=L,C,R 



+ V n , (1) 



where H a — ^(u a ) u a + ^(u a ) K a u a represents coupled 
harmonic oscillators, u a is a column vector consisting of all 
the displacement variables in region a, and u a is the corre- 
sponding conjugate momentum. K a is the spring constant 
matrix and V LC = (V CL ) T is the coupling matrix of the 
left lead to the central region; similarly for V . The dy- 
namic matrix of the full linear system is 



K 



K L yLC o 
yCL R C yCR 

V RC K R 



(2) 



There will be no interaction between the two leads. The 
nonlinear part of the interaction will take the form 

t r lY^rp C C C , 1 rp C C C C tn\ 

Vn = 3 1^ T ^ k U i U J U k + 4 T ^ kl U i U 3 U k U l ( 3 ) 
ijk ijkl 



for perturbation expansion, but it can be arbitrary. No- 
tice that the notation is valid independent of the actual 
physical dimensions used. 

For easy treatment, we consider the leads as quasi-one- 
dimensional lattices characterized by three square matri- 
ces of equal size, fcoo, fcn, and fcoi- The matrix fcoo is the 
block immediately adjacent to the center, while fcn, and 
koi = fcfo are repeated for the semi-infinite chain of the 
lead. They form the whole dynamic matrix of the lead in 
the following way: 



K R = 



fk 00 k 01 ■ • ■ ^ 

fcio fcn fcoi 
fcio fcn h i 

V • • • fcio " • • J 



(4) 



If there are some irregularities in the junction structure, 
we can enlarge the central part, so that the lead is always 
a regular lattice. The nearest neighbor nature of K is 
not a real limitation, since one can use a large supercell 
so that only the neighbor cells have interactions. If the 
cross sections of the leads are large, we can also model 
bulk walls. 

The semi-infinite nature of the leads is important con- 
ceptually. First, the heat bath by definition must be suffi- 
ciently large so that any finite energy transfer does not af- 
fect its temperature. Second, phonon waves scattered back 
into the bath will get dissipated and will not be reflected 
back to the central region. Any finite harmonic system will 
not have such a property. 



2 Ballistic thermal transport 
2.1 Landauer formula 

As early as 1957, Rolf Landauer presented a very intuitive 
interpretation of electron conduction in nanoscale junc- 
tions based on the concept of wave scattering [14|15j . The 
same argument can be applied to the phonon transport 
as well |16|17|18I19T20] . Heat current flowing from left to 
right through a junction connected to two leads at dif- 
ferent equilibrium heat-bath temperatures Tl and Tr is 
given by the Landauer formula 



/= r^.?uvT[uj](f L -h 

Jo 27r 



(5) 



where f L>R = {exp[hu/(k B T LtR )] - l} 1 is the Bose- 
Einstein (or Planck) distribution for phonons, and T[u>] is 
known as the transmission coefficient (or transmittance). 
The formula describes the situation that the central re- 
gion is small in comparison with the coherent length of 
the waves so that it is treated as purely elastic scattering 
without energy loss. The dissipation resides in the heat 
baths. We refer to such a situation as ballistic thermal 
transport. Note that the transmission coefficient is inde- 
pendent of the temperature. All the temperature depen- 
dences are in the distribution functions Jl and f R . 
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Fig. 1. (a) Dispersion relation Tiw vs. normalized wave number 
of a (6,0) carbon nanotube. (b) The phonon energy hu) vs. the 
transmission coefficient T\uj\. 



The thermal conductance is defined as the limit 



lim 

^T,T R - 



>T T L - T R 



(6) 



The connection with the usual definition of thermal con- 
ductivity in the Fourier's law, J — — kVT, is, k = aL/S, 
where L is the length of a system in the direction of heat 
flow, while S is the cross-section area of the system. For 
quasi-one-dimensional systems like carbon nanotubes, a 
is a better quantity to use, since the cross-section area S 
is not very well defined. 

One very good use of the Landauer formula is that it 
provides an upper bound for the thermal conductance of 
quasi-one-dimensional periodic systems |21j . This is ob- 
tained when all modes of waves are transmitted without 
scattering or reflecting back to the heat bath. In this ideal 
situation, the transmission coefficient consists of a set of 
steps as a function of uj taking only integer values, that is 



T[lu] — number of modes at frequency to. 



(7) 



Figure da) shows the dispersion relation of a carbon nan- 
otube of chirality (6,0) with the force constants computed 
from Gaussian 03 [22j . We optimize a 7-period, hydrogen 
terminated structure. The basis set is 6-31G, and the the- 
oretical model is b31yp DFT method. Once the dispersion 
relation is computed, the transmission coefficient can be 
obtained by counting the numbers of positive wavevector 
branches of phonon vibrations. Actually, the transmission 
curve in Fig. [IJb) is calculated by a more general NEGF 
method, which will be discussed in Sec. 13.21 

In the low-temperature limit, the high energy optical 
modes will not be occupied due to the Bosc-Einstein dis- 
tribution. Only the lowest energy, long wave modes are 
important. In this low-temperature limit, we have T[lo] « 
N m = 4 for any quasi-one-dimensional systems possess- 
ing translational invariance in three directions and rota- 



tional invariance along the axial direction (the four acous- 
tic modes correspond to one longitudinal mode, two trans- 
verse modes, and an extra 'twist mode'), and the thermal 
conductance is given by 



— hui N m -J- 
2tt dT 



N„ 



n 2 k 2 B T 



(8) 



independent of the material parameters. This is known as 
the universal quantum thermal conductance |17) . in anal- 
ogous to the universal electron conductance 2e 2 /h. The 
material constant independence stems from the fact that 
in one dimension, the product of phonon density of states 
and the group velocity is one. In higher dimensions, we 
must integrate in the fc-space, instead of just a single uj 

23J . The group velocity enters explicitly. The universality 
appears to have a deep connection to the maximum infor- 
mation or entropy flow in a channel [24125] . The approx- 
imation T[lu] w 4 is accurate when fc^T is much smaller 
than the lowest optical mode energy TiAui. For typical car- 
bon nanotubes, this is valid when temperature is of order 
of few kelvin to 10 K. For Si nanowires with transverse 
dimension of few hundreds nanometers, it is below I K. 
The spacing Aus can be estimated as Alo ~ c/y/S where c 
is sound velocity, and ^/~S is an estimate to the transverse 
dimension. Indeed, Schwab et al have done a remarkable 
experiment that verified the above theoretical prediction 

26J. This thermal conductance quantum is independent 
of the types of energy carriers; it is also the same for elec- 
trons [27], even photons |28J. 



2.2 Calculation of transmission coefficients 

When the transmission is not ideal, we must do a hard 
calculation for the transmission coefficient T[u>\. Several 
approaches are possible. Such considerations also justify 
the Landauer formula itself. Blencowe in ref. [19] presented 
a derivation using elastic wave models. In this article, our 
emphasis will be on atomistic models. We will give an out- 
line of a derivation using lattice models. Alternatively, the 
Landauer formula can be derived also from an NEGF for- 
mulism. But before doing this, we comment on the elastic 
models. 

Elastic wave continuum models have been used ex- 
tensively in connection with the universal thermal con- 
ductance |lb'll7ll8U9l2"0"] . This is justified on the ground 
that at very low temperatures, only the long-wave, low fre- 
quency acoustic phonons contribute to the thermal trans- 
port. Since the wavelengths are much larger than the lat- 
tice spacings of the materials, continuum models are good 
approximations. Often times, one takes a very simple scalar 
wave governed by the equation 



d 2 & 



c 2 V 2 <f = 0, 



(9) 



with appropriate (von Neumann) boundary conditions. 
The elastic wave theories have been used to investigate 
the robustness of the universal thermal conductance in 
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normalize the polarization vector such that e 
that 



various geometries including the effect of surface rough- 
ness [29130] , abrnpthinctions [31132] , bends [33] , and other 5 n _ n ,. The dagger stands for Hermitian conjugate. Note 
channel shapes [34135136] . Many papers appear recently on that e„, fc in general is complex. We demand that e* k = 
scalar elastic- wave models with various geometries such as e„ A particular vibrational mode is considered to be 
superlattices, T-junctions, stub junctions, etc |37|38|39|4U|41^)4^|M^5^|4g|^4^5p neral motion can be expressed 
It is important to ask the question, at what tempera- i n ' terms of normal mode coordinates Q n ,k(t), as u t (t) = 



ture do the elastic models break down? In an elastic model 
for a 2D channel, the dispersion relation takes the form 



u) 2 = (cqf + bn 2 



0,1,2, 



(10) 



where c is sound velocity, q is wave number, and b is some 
constant. Clearly the concept of Brillouin boundary is ab- 
sent in elastic systems. Also, the optical phonons are not 
represented correctly. Thus the elastic models will not be 
correct quantitatively when ksT > HAlo. When the typi- 
cal phonons have the wavelengths comparable to the lat- 
tice spacings of a crystal, we expect that the continuum 
theory breaks down. This happens when he /a ~ fcgT, 
where a is lattice constant. Alternatively, real crystals 
have a upper cut-off frequency, u) max , and associated De- 
bye temperature, Tjj — hu> max /kB- The values of temper- 
atures of the above estimates for real systems appear to 
be in the range 100 to 1000 K. It is then interesting to 
see a good comparison between elastic models and lattice 
models (e.g., carbon nanotubes) for thermal transport. 



2.2.1 Scattering picture and Landauer formula 

In this subsection, we give a derivation of the Landauer 
formula from an atomistic model of subsec. 11.11 following 
the steps of refs. |19|23|51|52j . First, in order to define the 
transmission amplitudes, we need to discuss the normal 
modes and the scattering problem. Let us consider only 
one lead as a quasi-one-dimensional system consisting of 
N repeating cells with periodic boundary conditions. We'll 
take N — > oo eventually. Due to periodicity, the displace- 
ments of the normal modes satisfy Bloch theorem, so that 
all the cells are oscillating in essentially the same way ex- 
cept a phase factor (see, e.g., ref. [53]). We write 



ui = eX e 



(11) 



where ui is a column vector for the displacements in cell I 
for M degrees of freedom, e is the polarization vector (of 
dimension M) satisfying the equation 



N 1 / 2 ^2 n k Q n ,k{t)ui } n,k- In terms of the normal modes, 
the total energy of the system is 



H 



R — n 



E 

n,k 



IQ„,fcWI 2 + < fc |Q„,fcWI : 



(13) 



In the derivation below, we need an expression for the 
group velocity v — duj/dq in terms of the eigenvectors. 
This is obtained by differentiating the eigenvalue equation, 
(112[) . with respect to the wavevector q, and multiplying the 
resulting equation with from left, 



2u>v 



1 



in e 



(14) 



We have used the Hermitian conjugate of Eq. (|12[) to elim- 
inate a term involving de/dq. 

In a scattering problem setting, we usually specify a 
frequency first, and ask the value of wavevector q for a 
given mode n. For this reason, we consider the normal 
modes as specified by 



Ul, n (w) = e„(w)e 



iq n (tu)la 



(15) 



We now introduce superscripts (a, a) to specify the side of 
the leads, a = L or R, and the direction of energy trans- 
mission a = ±. We'll use the convention that energy flow 
into (group velocity toward) the center is '+'. The left 
lead and right lead can be different lattices and thus have 
different dispersion relations and eigenmodes. For nota- 
tional simplicity, we'll omit the frequency argument, and 
use both (n, k) and (n, a, u>) to specify the modes inter- 
changeably. If a particular mode u\ L ' + ' is sent from the 
left lead to the central region, part of it will be reflected 
back, and part of it will be transmitted. For the cells far 
away from the center, the waves take the asymptotic form 



~(L, 

= U Ln 



on the left lead, and 



LL 



I 



(16) 



fcii-ifcio-AA dl )e = 0. (12) 



,RL 



(17) 



The identity matrix /, the spring constant matrices for a 
single cell, k\\, fcio, and fcoi, are size M x M. Imposing 
the Born- von Karman boundary condition, ui — w;+at, one 



gets A 



a i2-rzk/N 



q = 2nk/(aN); a is the lattice 



constant. We can choose the integer k — —N/2, —N/2 + 
1, • ■ • , N/2 — 1 for the linearly independent solutions. For 
each k, Eq. (JT3J) gives M eigenmodes with angular fre- 
quency uj n ^k and polarization vector where n = 1, 2, • • • 
M labels the branches, and k labels the wave vectors. We'll 



on the right lead. Equation (flTj) defines the transmission 
amplitude t^P n at frequency u> of the mode n to mode n! 
on the right, while Eq. (| 16|) defines the reflection ampli- 
tude back into mode n' . Similarly, we can send waves from 
the right lead to the left lead and define analogously t LR 
and t RR . From these transmission amplitudes, the total 
transmission, T[w], can be calculated. 

Our next task is to find a classical expression for the 
energy current and then to quantize the expression and 
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find quantum energy current. An expression for micro- 
scopic energy current can be obtained by considering the 
conservation of energy [9]. For a quasi-lD system, this is 

(18) 



dt 

where we define the local energy in cell I as 

H l = 2 (""f '"i + u T-ikoiUl + ujk lx ui + uf +1 k 10 ui) , 

(19) 

such that J2i Hi = Hr (or Hl). By differentiating Hi with 
respective to time t and using the equation of motion, we 
can see that 

h = ]^uj kxoui-x - ^uj^kmiii (20) 

satisfies the requirement. /; is the energy current from cell 
I — 1 to cell I. Expressing a general vibration as superpo- 
sition of modes with amplitudes Q n ,k, 



Ul(t) = /jr E Qn,k € n,k^ 
n,k 



+ C.C., 



(21) 



where c.c. stands for complex conjugate, and substituting 
it into Eq. (|20[) , and performing a time average, we obtain 



1 



1 



(II) = -^2^ iuj n,k\Qn,k\ e„, fe I Afcoi - jk 10 J e«,fc- (22) 

n,k ^ 

In deriving the above expression, we used the fact that 
the time average of e l ^" t: ~ w "' fc '^ is zero, unless n = n' 
and k — k' . The expression in the brackets can be further 
simplified in terms of the group velocity v n .k — duj n ^/dq. 
The final expression for the classical energy current in 
terms of the normal mode amplitudes is 



/ = 



—( 1 

2aN 

n,k 



J n,k 



\Qn,k 



(23) 



We now derive a set of constraint equations for the 
transmission amplitudes t"", due to energy conservation. 
Consider an arbitrary linear combination of excitations at 
frequency oj with incident waves from left and right, the 
total vibration at lead a can be written as 



£< 



a~(a, + ) 



"l.n 



E 



t aa ',a a ',v ( - a '~ 



(24) 



In the above formula, we have the convention that rep- 
resents the magnitude and thus is always positive. De- 
manding that II + Ir = for any amplitudes a", we 
obtain 

t f vt = v, (26) 

where t is a matrix with elements t"", where (a, n) are 
considered row index and (a', n') the column index, v is a 
diagonal matrix with the elements v" = i)"/o„ arranged 
in the same order as t. cll and an are the lattice constants 
of the left and right lead, respectively. The matrix t is 
somewhat close to be unitary, but is not. If we define S = 
v 1 / 2 tv -1 / 2 , then S is unitary. From SS^ = S^S = I, we 
can also show that 



tv~ 1 tt=v- 1 . 



(27) 



We now discuss the quantization of the problem. First, 
we consider only an isolated lead with periodic bound- 
ary conditions. Let us introduce the annihilation opera- 
tor in Heisenberg picture a rit k(t) — a n ^e~ lu>n - kt associ- 
ated with mode (n, k) and its Hermitian conjugate a k (t) 
for the creation operator, satisfying the usual commuta- 
tion relations, [a n> k, a„',fc'] = 0, [a^ fc , a) n , k ,] = 0, and 

[a n ^k, a} n , k ,] = Sn^n'Sk.k'- Then the canonical coordinate 
operator is 



Qn,k(t) = 



2ui 7} 



a„, k (t) + al k {t) 



(28) 



satisfying [Q„,fe(*), Ql',fc'(*)] = i^n,n'h,k' • Using the re- 
lation between the normal mode coordinates and the orig- 
inal coordinates, we can write 



hi 



(') = £■ 



£ n ,ke N a n ,k{t) + n.c, 



n , k 



(29) 



where h.c. stands for Hermitian conjugate of the proceed- 
ing term. We have [ui(t), uj,{t)] = ihlSij/, as required. 



We now substitute Eq. (|29|) into the current expression, 
Eq. (|20|) . taking a time average, and using the result of 
Eq. ([HD, we get 



1 



'' ' ,k 



(30) 



The first term above represents the superpositions of var- 
ious incident modes, while the second term contains re- 
flected waves, as well as incident waves from the opposite 
side. Comparing with the definition of normal mode ampli- 
tude, we find Q n ,k = VNa" for these n and k correspond- 
ing to incoming waves, and equal to \/N J2 n ' a' C a n' f° r 
these of outgoing waves. Substituting the values of Q n> k 
into Eq. (|2"3"|) , we find the incoming current from lead a to 
the center, 



I a (uj) 



2^(kMEc>"<'i 



(25) 



This is a well-known result for the energy current for a 
lattice [54] . 

For the open system with two leads and a central 
junction, we need to consider the two sides simultane- 
ously even if our interest is only the current on the left 
lead. Specifically the creation and annihilation operators 
will be associated with the scattering state Eq. (fl6)l and 
(fP7|) . originating from the left lead, a^ k , and a similar 

one originating from the right lead, a^ k . The general dis- 
placement operator u is a superposition of the scattering 
waves similar to Eq. ([24l . with the number a" replaced 
by the operator yjTi/ (2u>N)a < ^ k (t) and summed over all 
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modes, plus its h.c. After substituting the operator u into 
the current expression, Eq. (|20|). we obtain many terms, 
quadratic in a™ k or a" k . Now at this point, great simplifi- 
cation can be made if we evoke nonequilibrium thermody- 
namic average. We assume that the scattering states are 
distributed according to equilibrium distribution at tem- 
perature Tl = l/(fcB/3i) for the mode originating from the 
left lead, and Tr for those originating from the right lead. 
This causes an asymmetry for the left moving and right 
moving phonon waves, thus generates a thermal current. 
The thermal average (• • •) can be thought of as taking the 
trace of the creation/annihilation operators with respect 
to a density matrix given by p cx e~^ LHL ~^ aHrt , with the 
following results 



\ a n,k a n',k') — 0, 
( a n.k a n'k' ) = 0, 



i a n,k a n'k') ~ $n,n' &k,k' &u,a> f a, 



(31) 
(32) 
(33) 



where f a is the Bose-Einstein distribution function at the 
temperature T a . Using the above results, Eq. (|14[) . and 
with some algebra, we get, after cancelling the terms cor- 
responding to the zero-point motions, 



A = e iqa becomes an eigenvalue. Let us introduce a new 
vector £ the same dimension as e, by fc 10 e = AC, then 



u 2 I— fcn 







= A 



koi 
I 



(36) 



where / is an M x M identity matrix. Equation ([36]) sets 
up a so-called generalized eigenvalue problem. Since our 
eigenvalue problem has a dimension of 2M, we'll have upto 
2M eigenvalues for A. One can prove that the eigenvalues 
come in pair of inverse of each other, A + A_ = 1. Most 
of the eigenvalues represent evanescent modes, but these 
with |A| = 1 represent traveling waves. 

It is convenient to know the sign of the group velocity 
without actually calculating it. By adding a small imagi- 
nary part to ui, that is, replacing it by u) + if], then none 
of the eigenvalues A will have modulus exactly 1. Consid- 
ering r\ as a small perturbation, we find for the traveling 
waves [58] 

0+ 



= 1-7?- 



I] 



(37) 



That is, the forward moving waves with group velocity v > 
have |A| < 1. This result will be handy for constructing 
the surface Green's function of the lead later. 



ft.) = E 



TlUn.kV" l. 



k>0,n 



£ |0'| 2 /a'K, fe ) 



(34) 



Using the properties of the matrix t, Eq. (|26p and ([2"7]) . we 
can make the expression more symmetric with respect to 
the leads. We finally take the limit N — + oo by replacing 
J2k v n k/( a <*Na) • • • by J dLu/{2ir) ■ ■ ■. We obtain the Lan- 
dauer formula, Eq. ([5]), with the total transmission given 
by [52], 



T\uo] = 



v R 



n,n' 



ZHL\i RL I 
-r, \ L n'n\ 



(35) 



The above formula has a very simple physical interpreta- 
tion. The total transmission at a given frequency is a sum 
of all possible incoming channels n from the left and all 
outgoing channels n' on the right. Each term is simply 
the ratio of the energy current of the transmitted wave 
u[^, ^tn^ri t° the incident wave u\ L ^ + \ Since energy must 
be conserved, each term cannot exceed 1. 



2.2.2 Generalized eigenvalue problem 

The eigen modes can be found by Eq. (fTTj) and Eq. ([12"]) . 
In a standard dispersion relation calculation, one specifies 
a wave number q first, and then solves Eq. (fT2"]l to find 
various eigen frequencies. This is inconvenient for solving 
a scattering problem. Here we like to specify a frequency 
u> first, and find q and associated polarization vector e. To 
do this |55|56|57j . we rewrite Eq. (12]) in a form so that 



2.2.3 Scattering boundary condition method 

In ref. [52], Wang et al. proposed a scattering bound- 
ary condition method to compute the transmission co- 
efficient t",". The method is similar but not quite the 
same as the mode-matching methods 59 55 60 used in 
electronic transport. The idea is very simple. We write 
down the linear equations for the center, and impose the 
scattering boundary condition to solve it. Let us name 
the supercells consistently so that the left lead will have 
/ = •••, —2, —1, 0, the center consists of 1 to N, and right 
lead is N + 1, N+ 2, • • •. Then the equations that have not 
been dealt with by the leads are 



no M -i + (^ - ^oo) u o 



V LC u c = 0, (38) 



-V CL u Q + {lo 2 - K c )u c - V CR u N+l = 0, (39) 
-V RC u G + (lu 2 - k R )u N+1 - k R u N+2 = 0. (40) 

We have assumed that the interactions between the center 
and leads reach only the first layer of the leads immedi- 
ately adjacent to the center. In these equations u c is un- 
known, but u-i, uq, ujv+i, and ujv+2 are set to the values 
given by Eqs. |T6]) and fT7|. with t%P n and t R h n treated as 
unknowns. Since the form of Eqs. (fT6|) and (fT7]l is only 
asymptotic when the evanescent modes decay to zero, we 
should define the center part to include sections of the 
periodic leads. Alternatively, we should also include the 
evanescent modes in the summation of n' in the boundary 
conditions, but they should be excluded when calculating 
the total transmission. 

However, if the evanescent modes are omitted, we have 
a linear system with too many equations and few un- 
knowns. We solve the system numerically with a singular- 
value-decomposition method [61j , which gives a very good 
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approximation to the original problem when the evanes- 
cent modes can be approximated to zero. 



2.2.4 Mode-matching 



the transmission ampli- 



given incident mode Uq = e t | L , 
tude to n' is computed by solving Eqs. ([38]) . ([39)) , and 
(|40| . then obtained from the n' component of a vector 



RL 



UN+l J 



(46) 



Ando [59] , Khomyakov et al. [55160] , and also Ting et al. 
[62 63] gave an elegant method for computing the trans- 
mission amplitudes which we outline briefly here. The im- 
portant idea is to eliminate U-\ and ujv+2 in Eqs (|38p to 
(|40|) so that they form a closed set of equations that can 
be solved uniquely. One first solves the generalized eigen- 
value problem of Eq. (|36p . then writes a general solution as 
linear combination of all the modes (including the evanes- 
cent modes). We divide the modes according to |A|; these 
modes with |A| < 1 (forward or right moving waves) will 
be designated +, and those of |A| > 1 (backward moving 
waves) as — . We add a small positive imaginary number 
irj to lo so that none of the modes will have |A| = 1 exactly. 
Then we can write for the perfect periodic leads 



M 

a— ± n—1 



E+A'a + + E~A l a" 



(41) 



where we have written in a more compact form by intro- 
ducing matrices. E + = (e+, ej, • • •) is an M x M' matrix 
formed by the column vectors of e+; similarly for E~ . 
A± is an M' x M' diagonal matrix with the eigenvalues 
A„.±, while is a column vector for the amplitudes of 
the modes. We normalize the vectors, (e^)^e^ = 1- If fccu i s 
not singular, then M' = M, the number of degrees of free- 
dom in a cell. However, if it is singular, then M 1 < M . We 
exclude the trivial solutions = in Eq. (|4"Tj) . The solu- 



tion can be split into + and 
From Eq. (j4"Tj) . we see that 



component, ui 



"l+s 



?+ /l'+s a + 



E + A^ 
E+ A^iE+y 1 E+ A l + a^ 
F+(s)u+ 



(42) 



where F+(s) = E+A S + (E+)- 1 . If M' < M, the matrix in- 
verse should be interpreted as the pseudo-inverse [(E + )^ E^ 
(E + )^ . F~(s) is similarly defined. 

With the F operators defined for both the left and 
right lead, we can eliminate u_i and u^+2 using the equa- 
tions 



2.2.5 Scattering and transfer matrices 



Transfer matrix methods have been used for computing 
the transmission coefficients in quasi-one-dimensional atomic 
models [64 65 66 67 68 69 70J as well as continuum mod- 
els 37 36 . It is efficient in that the transfer matrix can be 
obtained in O(N) computer time where N is the size of 
the central region. However, if there are evanescent modes 
with large |A|, the evaluation of the transfer matrix can 
be numerically rather unstable, particularly if N is large. 
There is method to overcome this problem for special cases 

Let a \ be the incoming wave amplitude vector [c.f. 
Eq. (|4"Tj) ] from the left lead, and a^ be the incoming wave 
amplitude vector from the right lead, then the scattering 
matrix relates the incoming waves to the outgoing waves 



l - 



t LL t LR 
t RL t RR 



(47) 



which is the t matrix introduced in subsec. 12.2.11 The 
transfer matrix T (in the eigenmode representation), on 
the other hand, relates the amplitudes of the modes of left 
lead to the right lead: 



^11 T~\2 
T~2\ T~22 



(48) 



LR _ 



It can be checked that the transmission matrix t 
T 2 ~^ , or t RL = {T~ X )ii ■ So, knowing the transfer ma- 
trix one can deduce the transmission amplitudes. While 
the scattering matrix t is always square, the transfer ma- 
trix T can be rectangular. The matrix inverse should be 
interpreted as pseudo-inverse in such a case. 

We assume that the system can be divided into (prin- 
cipal) layers or cells, labelled by an integer index I, such 
that only the nearest neighbors have interactions. Then 

I the equations of motion at a fixed frequency lj are given 

'by 



kid-iui-i + (hj ~ lo 2 )ui + kij +1 u l+ i = 0. 



(49) 



7 l + (-l)<+^(-l)%> 



Uq = U(f + U Q , 

u N+2 = F+(l)u+ +1 = F+(l)u N+1 . 



(43) 
(44) 
(45) 



If the coupling matrix fc;,2+i is not singular, we can map 
from the displacements of ui-\,ui to ui, 



Ul+l 

ui 



Pi 



For the scattering state, there is no backward moving 



waves on the right, so u 



N+l 



0. We can eliminate 



u, 

Ul-l) 

.-1 , ,2 
"1,1 + 







Uj 
Ul- 



.(50) 



in favor of Uq and uq. With that, we move the terms in- 
volving Uq to the right-hand side of Eq. p8|) . For each 



We call Pi also the transfer matrix (or promotion matrix). 
The motions of the left lead and right lead can be now 
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related through 



2.2.6 NEGF: Caroli formula 



UN+2 

un+i 



= P(N+l : 0) 



U-i 



Another efficient way of computing the transmission coef- 
ficient is through the so-called Caroli formula: 



= P N+1 P N ---P 1 P Q 



u 



(51) 



tm = Tv(G r r L G a r R ), 



(57) 



We note that the transfer matrix for the leads, Pi when I < 
or I > N+l, becomes independent oil due to periodicity 
of the lattice; let us call them P. P has some interesting 
properties: det(P) = 1, and the eigenvalue problem of P 
is essentially the same as that of Eq. (f36|) . We can relate 
Eq. ([5"Tj) with (|48|) and obtain the transmission coefficient 
if we relate ui to the amplitude a. This is given by Eq. ([4"T]) ; 
rewritten in matrix form, we have 



m \ ( E+A l + E-A l _ 
ui-i) ~ \ E+A 1 - 1 E-A 1 - 1 



(52) 



Let us denote the transformation matrix from a to (ui, u;_i 
by Ul and Ur for the left and right leads, respectively, 
then we have T = f/^ x P(JV + 1 : 0)U L . On the other 
hand, we could also use Eq. fl3|. (|44|) . and (gSj) to elimi- 
nate ujv+2, and solve Eq. (|5ip for ujv+i and obtain 
the transmission coefficient via Eq. (|46|) . 

An interesting relation exists between the eigenvalues 
of 'square' of the transfer matrix and the total transmis- 
sion due to Pichard [72]. The results are first derived for 
the electronic transport where the scattering matrix is uni- 
tary. Since in phonon transport t is not unitary, we cannot 
apply the formula directly. Instead, we need to work with 
S = v^tv -1 / 2 , and replace the amplitude vector a by 
c = v 1 / 2 a. Then we have 



S 



,+ 



(53) 



"R/ \YR, 

The corresponding transfer matrix associated with S is 



v 



~l/2 




r,-V2 



v 





r -i/2 



(54) 



With S and M we have the Piccard identity [73174] 



[21 + MM ] + {MM^Y 1 } =- 



1 (SblSU 



4 V 







(55) 

where S aa > = Va t aa v a ?~ , a — L,R. The eigenvalues of 
the product AiAi^ are real and positive and come in pair 
of inverse of each other. We denote them by exp(±2xfc). 
Taking the trace of both sides of Eq. (l55|) , and using the 
fact that the total transmission T[u>\ = Ti[SRLSn L ), we 
obtain 

r N = E^- (56) 

, lush dj 

In applying the above formula, we should use a global 
transfer matrix P(N + 1 + s : — s) for a sufficiently large s 
to "damp out" the contributions from evanescent modes. 



where G r = (G a )' f is the retarded Green's function for 
the central region, while Pl,r describes the interaction 
between the leads and the central region. Caroli et al. 
first obtained a formula for the electronic transport in a 
slightly more restricted case [75] . Meir and Wingreen gave 
the above form from NEGF formulism [76]. For thermal 
transport, it has been derived by many authors from var- 
ious perspectives [771781791801811821831121131 . 

In this subsection, we only discuss the relation of Eq. f5"7l 
with the surface Green's functions. We'll consider a deriva- 
tion in sec. 13.21 The retarded Green's function of the full 
linear system, including the leads, in frequency domain is 
given by the solution of the equation [(w + irf) 2 I — K\ G = 
yrl. Both K and G are (infinite) matrices indexed by the 
labels of atomic degrees of freedom. Partitioning the ma- 
trices according to the left, center, and right regions, as 
given in Eq. @ for K, and similarly for G, we have [32 84 

\uj+irffl-K L ~V LC 
-V CL (LU+ir,) 2 I-K c -V CR 
-V RC (u;+ir 1 ) 2 I-K R J 

(G LL G LC G LR \ //0 0\ 
x ( G CL G cc G CR ) = [ I ) . (58) 
\G RL G RC G RR j \00IJ 

Considering only the three equations formed by each row 
of the first matrix multiplying the middle column of G 
matrix, we have 

[(w + ivfl ~ K L ] G LC - V LC G CC = 0,(59) 

_yCL G LC + [( w + lI) )2 / _ ]f C] G CC _ yCRQRC = j {6Q) 
_yRC G CC + [ (w + ir)) 2j _ R R ]G RC = (61) 

Introducing the retarded surface Green's functions for the 
leads 

g r a =[{u + irifl-K a Y\ a — L,R, (62) 

and using Eqs. flS5J| and ([HI]) to eliminate G LC and G RC , 
we can express the central part of the Green's function as 

G r = G cc = [(w + i n fl ~K C -Sl- E R ] -\ (63) 

where retarded self-energy of the leads is given by 

K = V Ca g r a V aC . (64) 

G r is the matrix entering into the Caroli formula. The r a 
functions are given by 

r a = i(U r a ~U^ = -2lmV Ca g r a V aC . (65) 

As we can see, the transmission coefficient can be obtained 
if we know the surface Green's function g r a . We'll discuss 
the algorithms for computing the surface Green's func- 
tions in the next section. 
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2.2.7 Algorithms for surface Green's functions 

In an NEGF calculation for the total transmission, an im- 
portant piece of code is for the computation of the Green's 
function of an isolated, semi- infinite lead, satisfying, e.g., 
for the right lead, 



[{u + in?i-K R ] g r R = i, 



(66) 



where K is given by Eq. Q. In block matrix form, for 
the left-most column of g r R , we have 

[(u> + iri) 2 I - fcoojffoo - fcoiffio = I, (67) 
-kiogi-i,o + [(w + ir/) 2 I - hi] 9l,o 

-koigi + i,o = 0, Z = 1,2, • • • (68) 

In order to be consistent with the literature on surface 
Green's functions, we label the right lead starting from 
layer I = (instead of N + 1 used in Sec. 12.2.3ft . Since the 
coupling from the center to the lead only reaches the 0-th 
layer, we only need the M x M block matrix goo- Various 
methods have been used |85|86|87|88|89|90|91j to find 500, 
see [55] for an recent review. 

An equation can be obtained for 300 by considering 
the relation of the original system K R and a new system 
with the 0-th layer removed. Similar to the derivation of 
Eq. ([53]). this gives 



.9oo = [(w + irj) 2 I - fc o - k 01 gki 



(69) 



where g is the surface Green's function for a system with 
the 0-th layer deleted, i.e., g is the top-left element satis- 
fying a similar equation as Eq. (|66p . but with a K matrix 
obtained from K R with the 0-th block row and 0-th block 
column deleted. If we continue this operation of removing 
the next layer, we see that the K matrix will not change 
(since it is semi- infinite and translationally invariant), we 
must have 



g = [(uj + irj) 2 I - fcn - k 01 gk 
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(70) 



The above equation can be iterated starting from g = 0. 
This is known as simple iteration method |58j . However, 
such iterations converge slowly. After i iterations, the re- 
sult is equivalent to a truncated system of finite length of 
i layers. 

A more efficient and popular method is that of Lopez 
Sancho et al. [88] , which doubles the number of layers in 
each iteration, i.e., after i iterations, the effect of 2 l layers 
is taken into account. We refer to ref. [88] for its derivation. 
Here we give a pseudo-code for such an algorithm. The 
left arrow ' below denotes assignment, e is an error 
tolerance, and | • | stands for some matrix norm. 

s <- fcoo 
e <- fcn 
a <— fc i 
do { 



0- 



[(u + ill) 2 1 - e] 

T 

a 



s «— s + agf3 
e «— e + ag(3 + 
a <— ago. 
} while (|a| > e) 
500 *- [(w + irf) 2 ! 



/3ga 



The most accurate and perhaps also the fastest meth- 
ods 85 89 90 91J are based on the solution of the gener- 
alized eigenvalue problem, Eq. (|36p . We notice that gi t o 
and the general displacement ui satisfy exactly the same 
equation, for / > 0. Thus, we can express the solution of 
gt t Q in the general form of Eq. ([4T]) . However, since g^o 
must be finite in the limit I — > +00, we can only have the 
decay modes (+ components) with |A| < 1. In particular, 
Eq. ([42]) is also valid for gi t o, that is, 



g l+sfl = F+(s)g li0 , F+(s) = E+AUE+y 



(71) 



Substituting this result for / = 
can solve for goo explicitly as 



0, s = 1 into Eq. flfffl), we 



ffoo 



[{LO + ir,) 2 I~k m -k 01 F+{l)\ 



(72) 



The eigenvalue problem and matrix inverses all have 0(M 3 ) 
in computational complexity. We note that when con- 
structing the F + matrix, we need to include all the |A| < 1 
modes, including these with A = when fcoi is singular. 
The surface Green's function of the left lead can be com- 
puted by the same routine by first inverting the matrix 
indices, applying the algorithm, and then inverting goo 
back again. 



2.2.8 Relation between scattering picture and NEGF 

The results for the transmission coefficient from the scat- 
tering picture, Eq. ([33]) . and the Caroli formula, Eq. (T57]) . 
look quite different. It is shown explicitly in ref. [60] that 
they are completely equivalent. This can be done by first 
solving formally the mode-matching problem discussed in 
sec. 12.2.41 in terms of Green's function, given 

un+i — G]sr+ifiQo£ niL , (73) 
Qo = *&[*£(-!)- *!(-!)]■ (74) 

The Green's function in the above formula is the top-right 
corner block of G RL in Eq. ([58]) . which is related to the 
central part Green's function G r by 



Gn+i,o - g T Rfio vR gtvCL 9l,oo- 



(75) 



We can relate Qo with the surface Green's functions of 
the leads, and relate group velocities to self-energies, and 
eventually arrive at [BU] 



t RL n = %2uj 



Gn+i,o {(e l y 



(76) 



We can think of this result as the Fisher-Lee relation [92 
for the ballistic phonon transport. 

To illustrate the methods, we consider a defect model 
of carbon nanotube where the mass of one of the carbon 
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Fig. 2. Partial transmission T n i n [u] of the TW (twist) mode 
for a (8,0) carbon nanotube with a pinning defect site. The 
curve with label TW is from TW to TW, and that with TA is 
from TW to TA, other unlablled curves correspond to TW to 
optical modes. The values given are sum of the two degenerate 
modes of TA or optical modes. 

atom is set to infinite, creating a strong pinning at the 
site. The force constants are obtained from Brenner po- 
tential [93] ; starting with a perfect nanotube structure, the 
values of rows and columns of the K matrix correspond- 
ing to the defect atom are set to 0. Unlike some of the 
other defects 181I94I95| . such pinning has a strong effect 
to the low frequency transmissions due to broken transla- 
tional invariance. It also produces strong mode mixings. 
In Fig. [5J we show the partial transmission, defined as, 
T„/„[w] = |^4l 2 ^n'/^n- This quantity can be computed 
by the mode-matching method of sec. l2.2.41 or by Eq. (fTrll) . 
We observe that substantial fraction of the twist (TW) 
mode vibration is converted into transverse acoustic (TA) 
modes and other optical modes. Due to symmetric reason, 
the longitudinal mode is not excited and has zero trans- 
mission. The peaks in the TW to TW transmission are 
associated with the presence of an optical mode with zero 
group velocity. 



2.2.9 Transmission from molecular dynamics 

Schelling et al. proposed a MD method to compute the 
transmission coefficient [96 97 98 99 94[. The idea is to 
send a wavepacket localized both in real space and mo- 
mentum space, and to measure the fraction of energy 
transmitted from one lead to the other. The method is 
applied to semi-conductor interfaces and nanowires, grain 
boundaries, as well as carbon nanotubes with point de- 
fects. If the amplitudes of the waves are small enough, the 
results are effectively the same as that of lattice dynam- 
ics or NEGF. Due to limited resolution in both space and 
time in MD, the results are typically less accurate. How- 
ever, the simulation can give a good visual presentation 
of the phonon dynamics. 



The initial wavepacket can be constructed either in 
real space or in momentum space. For example, we can 
consider the superposition of eigenmodes of different fre- 
quencies, 

/ + OC 7 
a n ^)el L (w)X l -i°^)e-^— + c.c. (77) 

If we choose a n {ui) cx e^""" - 1 \ we'll have a wave- 
packet centered around Zo at t = in real space with a 
spread of order v n /a lattice spacings. 



2.3 Dimension crossover, etc. 

In this section, we mention some of the work that has 
been done for ballistic phonon transport, using the tech- 
niques mentioned earlier. One interesting aspect of ballis- 
tic transport is the dimension effect. For quasi-one-dimen- 
sional systems, at sufficiently low temperatures, we expect 
that the conductance depends on temperature linearly, 
a cx T . Calculations on carbon nanotubes and graphene 
strips [1001101] confirm this expectation. It is interesting 
to know, but to our knowledge not answered, how the 
ID behavior crosses over to a full 2D behavior. Naively, 
we expect that a cx T d , where d is the dimensionality. 
In two dimensions, this is not the case, but it is given 
by T 3 / 2 [21|102j . This is because that one of the acous- 
tic branch has a dispersion relation that is quartic in the 
wavevector. In ref. |103j . diameter dependence of silicon 
nanowires is studied using NEGF method. It is found that, 
as the diameter increases, the effective exponent increases 
from 1 to 3, where d = 3 being the bulk value. Interfaces 
[104] and atomistic junction systems [105ll06j have been 
studied. The effect to thermal current of Stone- Wales de- 
fects in graphene nanoribbons is studied in [107] ; and the 
phonon transmission through defects in carbon nanotubes 
using first principles force constants is studied in [108] . 
In ref. [109j Landauer-like formula for thermal transport 
through a mesoscopic weak link is given. 



3 Nonlinear effects in thermal transport 

3.1 Phenomenological treatment of nonlinear 
interactions: effective transmission 

If we take a short piece of quasi-one-dimensional material, 
such as a carbon nanotube, with an ideal lead-system con- 
tact and ignoring nonlinearity, the thermal energy trans- 
mission is given by the Landauer formula with the trans- 
mission coefficient given by Eq. (|7|). Each mode n is de- 
scribed by a simple rectangular function, T„[o;] = 1 for the 
frequency lu within the phonon band of the n-th branch, 
and T„[w] = outside the band. However, this picture 
breaks down if the system is long in comparison with the 
phonon mean free path. In ref. [HOj (and also [70]), a phe- 
nomenological description of the effect of finite phonon 
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Fig. 3. Carbon nanotube thermal transport in ballistic- 
diffusive regimes: (a) Length dependence of thermal conduc- 
tivity for different chirality SWCNTs at 300 K; (b) Average 
mean free path for nanotube (6,6) as a function of tempera- 
ture; (c) Change of the effective power-law exponents for short 
and long CNTs with temperature; (d) Change of the exponents 
with diameter. From ref. I110I. 



mean free path is taken into account in analogy to the 
electron transport [lllj : 



T n [u] = (l + L/£ n (u)Y 



(78) 



where £ n (uj) is the mean free path of phonons of mode n 
with frequency uj (A quite different form, l/Tfo/j oc (u — 
loq) 2 + d 2 , was used in ref. [112] . which cannot account for 
the nonlinear effect). This formula interpolates between 
the ballistic transport and diffusive transport. When l n ^> 
L we obtain the Landauer formula, while in the opposite 
limit, we reproduce the well-known Debye-Peierls formula 
for the thermal conductivity. 

The mean free path is related to the phonon relaxation 
time by £ n = r n v n , where v n is phonon group velocity. The 
relaxation time can be obtained by considering the de- 
tailed scattering mechanisms, such as Umklapp processes, 
defect and boundary scatterings [1131114111511 1611 1711 18IT 
When detailed structure and nonlinear information are 
put in, Eq. 1(75]) should be able to give quantitatively cor- 
rect predictions. In ref. [110] . a phenomenological expres- 
sion I w A/(oj 2 T) is used for a calculation of carbon nan- 
otube. The divergence thermal conductivity exponents are 
defined by k oc L as or L° L for short and long carbon nan- 
otubes. As shown in Fig. [31 at a fixed temperature, the 
thermal conductivity grows with system sizes with a large 
exponent for short nanotubes and reduces to a smaller 
exponent for much longer tubes. We notice that purely 
diffusive behavior is not observed even for systems as long 
as 100 jim. 



3.2 Nonequilibrium Green's function method 

The nonequilibrium Green's function (NEGF) method is 
an elegant and powerful method to treat nonequilibrium 
and interacting systems in a rigorous way. Its root is from 
quantum field theory. Some of the early formulations are 
by Schwinger [121] , Kadanoff and Baym [122] . and per- 
haps most significantly by Keldysh [123] , The books of 
refs. [1241125] contain some interesting historical aspects 
of the developments. The Keldysh diagrammatic expan- 
sion method is also generalized to cases of arbitrary (cor- 
related) initial states [126] . Meir et al. [761127] applied 
the technique to electronic transport through junctions. 
We recommend the books by Datta [111] and Haug and 
Jauho [128] for some necessary backgrounds in this topic. 
The application of NEGF method to thermal transport is 
relatively recent. In refs. [77I129I8T] . ballistic transport is 
treated. However, the real strength of the NEGF method 
is its ability to handle nonlinear interactions rigorously, at 
least in principle [12113113011311531 . 



3.2.1 Definitions of the Green's functions and their relations 

In NEGF approach, it is convenient and also necessary 
to introduce several different versions of the Green's func- 
tions. We will start with the definition of six Green's func- 
tions [1321133] : 



G r (t,t') = -i6(t-t'){[u(t),u(t') T ]), 

G>(t,t') = -i(u(t)u(t') T ), 
G<(t,t') = ~i(u{t')u(t) T ) T , 



(79) 
(80) 
(81) 
(82) 



G*(i, t') =6(t-t')G > (t,t')+6(t' -t)G < {t,t'), (83) 
G'(M') = 0(t' - t)G> (t,t') + 6(t - t')G<(t,t'). (84) 

They are known as retarded, advanced, greater, lesser, 
time-ordered, and anti-time ordered Green's functions, re- 
spectively. u(t) is a column vector of the particle displace- 
ment in Heisenberg picture. The step function 6{t) = 1 if 
t > and if t < 0. The notation ([A, B T ]) represents 
a matrix and should be interpreted as (AB T ) — (BA T ) T . 
To avoid the Ti floating around, we have set Ti — 1 in this 
section. 

In equilibrium or nonequilibrium steady states, the 
Green's functions depend only on the difference in time, 
101128] . The Fourier transform of G r (t - t') = G r (t 7 t') is 
defined as 



G r [u>] 



G r {t)e iujt dt. 



(85) 



Note that we use square brackets to delimit the argument 
for the Green's functions in frequency domain. The fol- 
lowing linear relations hold in both frequency and time 
domains from the basic definitions: 



G r -G a = G> -G<, 
G t +G i = G > +G < , 
G 1 - G l = G r + G a . 



(86) 
(87) 
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Also, the relations G r = G* - G< and G a = G< - G l 
are useful. Out of the six Green's functions, only three of 
them are linearly independent. However, in systems with 
time translational invariance, the functions G r and G a are 
Hermitian conjugate of one other: 



H,+H r +H„+V +H„ 



G>] = (G>])t. 



(89) 



So in general nonequilibrium steady-state situations, only 
two of them are independent. We consider them to be G r 
and G < , but other choices are possible. There are other 
relations in the frequency domain as well: 

G < M t = -G < H, (90) 
G r [-u] = G r [uj]\ (91) 
G< [- u] = G> [lo] t = -G< [lo] *+ G r [uf- G r [lo]*. (92) 

The last two equations show that we only need to compute 
the positive frequency part of the functions. 

Equations (f56"|) to (f9"2")) are generally valid for nonequi- 
librium steady states. In thermal equilibrium, there is an 
additional equation relating G r and G < : 



G<H = f(co) G>]-G>] 



(93) 



where f(w) is the Bose-Einstein distribution function at 
temperature T. Equation |93|) is obtained by writing the 
Green's functions as a sum of energy eigenstates (known 
as the Lehmann representation, see, e.g., ref. [134] ). In 
equilibrium, we also have G > [w] = e^ 3ul G < [uj]. Thus in 
equilibrium, there is only one independent Green's func- 
tion; we take it to be G r . 

The last definition of Green's functions is the contour- 
ordered Green's function, which is the most convenient 
entity for a diagrammatic expansion. It is defined by 



G(t,t') = ~i(T t u(t)u(t') t ), 



(94) 



where the variable r is on a Keldysh contour from — oo 
to +oo and back from +oo to — oo. The contour-ordered 
Green's function includes four different Green's functions 
given earlier: 



G T<T (t,t') = lim G(t + iea,t' + iea'), 



a = ±(l). (95) 



We have introduced a branch index a, such that r = t + 
iea. This is purely a book-keeping device to indicate the 
branch that the r variable is at, a = +1 means r is at 
the — oo to +oo branch, while a = — 1 means r is at the 
returning branch. With this notation, we can identify that 



C++ = G', G~ 
in matrix form 



G\ G+- = G< and G~+ = G> 



G(r,r') 



G* G< 
G> G* 



(96) 



In dealing with the contour-ordered Green's functions, 
we often encounter convolution of the form 

B{t,t') =JdT 1 Jdr 2 ■ ■A 1 (t,t 1 )A 2 (ti,t 2 ) ■ •A n (T n -i,r'). 

(97) 




Equilibrium at T a 



t = 
Nonequilibrium 
steady state 



Fig. 4. A schematic to illustrate the two adiabatic switch-ons. 
From ref. Q3]. 



This form of expression can be easily translated into the 
retarded and lesser Green's functions in frequency domain 
by the Langreth theorem as |135|128|136j 



(98) 



(99) 



B r [u>]=A\[w)Al[w)---A r M, n = 2,3,... 
B<[u}=Al[u;}---A r n _ 1 [uj}A<[u;} + 

A r M---A r n _ 2 [u]A<_MA a M + 
3.2.2 Green's functions of the free leads 



We will adopt the following notations, the Green's func- 
tions of the semi-infinite, equilibrium free left or right lead 
will be denoted <?, while the ballistic system of nonequi- 
librium, noninteracting Green's functions will be denoted 
Go, and the full interacting system will be denoted by G. 
Their relations are shown schematically in Fig. 2] through 
two adiabatic switch-on's of the interactions. In this sub- 
section, we discuss the Green's function of the leads. To 
make a connection with the Green's functions discussed 
in sec. 12.2.61 we consider the equations of motion of the 
Green's functions of the right lead. By differentiating the 
definitions of various Green's functions with respect to 
time t twice, using the Heisenberg equation of motion of 
any dynamic variable A(t) of the right lead, 



i^- = [A(t),H R ], 
we can obtain the following set of equations, 



(100) 



d 2 g r '°" t (t,t') 

dt 2 



+ K H g™*(t, = -IS(t - t'), (101) 



at 2 

d 2 g < ' > (t,t') 
dt 2 



+ K H g t (t,t') = IS(t-t'), (102) 
+ K R g <> {t,t') =0. (103) 



We can also give an equation for the contour ordered ver- 



Ot 2 



+ K g(r, t') = -IS(t,t'), 



(104) 



where the r variable takes values on the Keldysh contour. 
The contour ordered version is related to the ordinary 
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Green's functions through Eq. |96|) . In order to be con- 
sistent with the equations for g and g', we must define a 
generalized <5-function on the contour as 



5{t, t 1 )^ a5^ a ,8{t-t'). 



(105) 



Solutions to these differential equations can be ob- 
tained by Fourier transform. For example, we have 



(u*-K*)g™*[u>] = L 



(106) 



Although g r , g a , and g l satisfy the same equation, their 
solutions are clearly different. They have different bound- 
ary conditions. The correct choice for the retarded Green's 
function g r is 



g r [uj] = [(uj + irj) 2 



K 1 



(107) 



where ry is an infinitesimal positive quantity to single out 
the correct path around the poles when performing an in- 
verse Fourier transform, such that g r (t) = for t < 0. 
Other Green's functions can be obtained through the gen- 
eral relations among the Green's functions, e.g., 3 < [w] = 

/MGfM -$>])■ 

3.2.3 Green's function Go and G 

To compute the Green's functions of the noncquilibrium 
and interacting systems, we need to use perturbation the- 
ory and the concept of adiabatic switch-on. We can think 
of two adiabatic switch-on's as illustrated in Fig. HI We 
imagine that at t = — oo the system has three decoupled 
regions, each at separate temperatures, Tl, To, and Tr. 
The couplings between the regions and the nonlinear inter- 
actions are turned off. The equilibrium Green's functions 
g a at temperature T a are known and take the form of 
Eq. (flQ7l) . The couplings V LC and V CR are then turned 
on slowly, and a steady state of the linear system is es- 
tablished at some time to <C 0. For this linear problem, 
the result does not depend on Tq', the initial condition of 
the finite center part is forgotten. Finally, the nonlinear 
interaction V n is turned on, and at time t = 0, a nonequi- 
librium steady state is established. 

The density matrices at time t = — oo, to, and t = are 
related in the following way in the (respective) interaction 
pictures: 

p(to) = S Q (t ,~oo)p(~ao)So(~oo,t ), (108) 
S (t,0 =Te- i ti v{t " )dt " , (109) 



p(0) = S(0,t o ) P (t o )S(t o ,0), 
S{t,t') = Te- i ti v ~VW, 



(110) 
(111) 



where T is the time-order operator (assuming t > t'). 

The contour-ordered Green's function can be obtained 
by a perturbation expansion of the interaction picture evo- 
lution operators (scattering matrix operators) and is ex- 
pressed as: 

G (r, t') = -t(T T u(T)u T (r')e- l I V{T " )dT ") g , (112) 
G(r,r') = -i{T T u{T)u T {r')e- i i V ^ dT '') Go ,{in) 



where in Eq. (]112p the unperturbed system is the uncou- 
pled system with Green's function g, while in Eq. (|113|) we 
make perturbative expansion from Go- Since the coupling 
V is quadratic, its expansion leads to exact results, e.g., 
Dyson equation for the central part of variables (dropping 
the superscript cc), 

G (T,T')=g C (T,T')+ fdT 1 dT 2 g C (T,T 1 )S(T 1 ,T 2 )Go(T 2 ,T'), 



where the self-energy due to linear interactions with the 
leads is 



S(t 1 ,t 2 )=V cl 9 l (t 1 ,t 2 )V 



tLC 



V 



CRR. 



(n,T 2 )V 



RC 



(H5) 

Thus, the Green's function Go can be expressed in terms 
of g exactly. Using the Langreth theorem discussed at the 
end of Sec. 13. 2.11 the contour ordered Dyson equation gives 
two independent equations, the retarded version has solu- 
tion identical to Eq. (]63p . while the lesser component can 
be solved to given [128) 



G<[u]=G r MZ<[Lo]G a M 



(116) 



However, for the nonlinear problem of the second ex- 
pansion, we need to use the machinery of Feynman dia- 
grammatic technique [13711381139] . Since the nonlinear in- 
teractions depend on the displacement u , we only need to 
consider the Green's functions G cc . By expanding the ex- 
ponential in Eq. (j!13p . a series in the nonlinear interaction 
strength is obtained. The Feynman diagrams generated 
here are identical to the field theories of (f> 3 and <fi 4 [140] . 
and are similar to that in refs. [14111421143] for phonon 
retarded self-energies. Since the structure of the diagram- 
matic expansion of the contour ordered Green's functions 
is identical to the standard ground state (T = 0) expan- 
sion for the time-ordered Green's function, the Wick's the- 
orem [138] is also applicable here. The result of each term 
in the expansion can be expressed as product of Go- 

The expansion contains connected as well as discon- 
nected Feynman diagrams. The disconnected diagrams are 
constant in time (not necessarily zero), and give rise to 
a thermal expansion effect. We can show that these dia- 
grams do not contribute to the thermal transport, as they 
are proportional to 5(u>) in the frequency domain. The 
thermal current formula has a factor of lu which makes it 
zero. The connected part of the Green's function satisfies 
a similar contour-ordered Dyson equation 133 relating 
G c to Go through a nonlinear self-energy U n : 

G c (r,r')=G (r,r')+JdT 1 dT2G (T,n)S n {n,r 2 )G c (T2,T'). 

In ordinary Green's functions and in frequency domain 
(to argument suppressed), the above Dyson equation has 
solutions 1281: 



Gl = ({u + l r 1 fl-K c -S^ K) 



(118) 
(119) 



The Eq. (j!19[) is known as the Keldysh equation. The 
connected part of the Green's function G c is needed to 
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f 1 + (-6) 



<3 *»qp 



Fig. 5. The Feynman diagrams for nonlinear self-energy E n 
with the prefactors for graphs of cubic and quartic interactions 
to second order in fi. From ref. 1131. 



compute the heat current. The extra subscript c will be 
dropped from this point. 



3.2.4 Equation of motion method 

An alternative method to obtain the Green's functions of 
interacting systems is through the equation of motion of 
the Green's functions. The equations of motion for various 
Green's functions may be used as a starting point to form 
mean- field type approximations. In ref. |144j the equation 
of motion for the nonequilibrium Green's function is jus- 
tified from the Keldysh formalism, so the two approaches 
are completely equivalent. 

It is convenient to work directly on the contour ordered 
Green's functions. To represent the ordering along the 
Keldysh contour, we introduce a generalized ^-function 
[145 , 6>(t, r'), which is defined to be 1 if r is later than 
t' along the contour and otherwise. Its derivative along 
the contour is the generalized ^-function, 



S(t,t') = 



86 (r, t') 
dr ' 



(120) 



The generalized (^-function is related to its (a, a') com- 
ponents by Eq. (|105[) . The retarded, advanced, and time- 
ordered components of S(r, r') is the same as the ordinary 
5(t—t'), lesser and greater component is 0, while anti-time- 
ordered component is —S(t—t'). We also need to generalize 
the canonical commutation relation on the contour, as 



i 6 



(121) 



Similarly, we also generalize the Heisenberg equation of 
motion, 

dA(r) 



ch 



= %[H(t),A(t) 



(122) 



Equations (|12ip and (|122[) are the same as the usual canon- 
ical commutation relation and the Heisenberg equation of 



motion on the upper branch (— oo to oo) of the Keldysh 
contour; but on the lower return branch, they are the time- 
reversed version of the equations (since dr = adt). 

We first consider the linear case, V n = 0. If we view the 
system as a whole, the contour ordered Green's function 
satisfies 



d 2 G (T,T>) 
dr 2 



-KG (t,t')=I5(t,t'). 



(123) 



This is obtained from taking derivatives twice to the def- 
inition of the contour-ordered Green's function, and us- 
ing Eqs. (|120p and (|12ip . If we partition the matrix Go 

to the submatrices Gq'" , a,a' — L, C, R, and similarly 
for K, we can obtain an equation similar in structure to 
Eq. (|58|) with uo 2 replaced by —I-§pz and I by I5(t, t') on 
the right-hand side of Eq. (|58p . Using a similar solution 
technique, we can solve for the contour-ordered nonequi- 
librium Green's functions 



G^(r,r') = / dr"g^ry)V LC G^{T"y) 1 (124) 



G^(r,r')=g c (r,r') + 



dn / dr 2 g c {t.t^S^^G^ {t 2i t'). (125) 



The first of the above equations relates the central part of 
the Green's functions to the lead-center Green's function, 
which, it turns out, is also valid when the nonlinear inter- 
action in the center is nonzero. The second equation is the 
Dyson equation relating the free leads to coupled linear 
system. The self-energy Z , (ti,t 2 ) is given by Eq. (|115l) . 

For the expansion of the Green's function G in terms 
of Go , we need to define a general n-point contour-ordered 
Green's function as 



G 



'(ti,t 2 ,---,t„) = 



.<»(r„)).(126) 



The index a = L,C,R labels the region, j labels the de- 
grees of freedom in that region, and r is the contour vari- 
able. The function is symmetric with respect to simulta- 
neous permutations of the triplet (a,j,r). 

We can write the interactions in a more symmetric 
form, e.g., for the cubic term, 



V n {r)dT - 



i j k 



T %3k {r, t',t") 
Up (r)uf (r')ufe (r") drdr' dr" . (127) 



Since the contour integral is J dr = a J_°° dt, we must 
define 

%(r,r',r") = T ijk 6(r, t')6(t, t") 

-> T m a'8 a<a/ S(t - t')a"5^ a „6(t - t").(128) 
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n n 
' ll ' l 

3/ S3 




Fig. 6. Recursive expansion rule for Green's functions. A ver- 
tex with n double lines denotes an n-point Green's function. 
The single line denotes Go- A single-line three-terminal vertex 
is associated with T(t,t' , t"). From ref. |13| . 



We now consider the equation of motion of the central 
part Green's function G cc when there is a cubic nonlinear 
interaction. Differentiating the Green's function with re- 
spect to r twice, using the Heisenberg equation of motion 
(on the contour r), we get 



Q T ,2 



+ Y. G ^{r,r')Kg = -5 jl 5{ry) 



tRC 
'kl 



k k 

£ // Gf C (T,T ll T 2 )r ill (T 1 ,r 2 ,r') ( ir 1( lr, (129) 



The Green's functions G CL and G can be eliminated 
in favor of G cc through (transpose of) Eq. (|124|) . Then 
the differential-integral equation for G cc can be solved in 
terms of G ccc and G% G with the help of Eq. ljI25)l . We 
obtain 

GT{t,t') = G^(r,r') + E /// G™(r,n) 



ikrn 



Tikm(n,T2,T3)G kml (T2,T3,T')dndT 2 dT3. (130) 

This equation can be represented by a diagram. Since G 
is related to a three-point Green's function, we also need 
the equation of motion of G ccc . We can derive general 
rules [13j for n-point Green's functions, summarized as 
follows (as shown in Fig. [6J: 

(1) Replace leg 1 by inserting a nonlinear coupling 
T(r, r' ,t") such that the outer leg is immediately con- 
nected with Go while the other two terminals increase the 
order of the Green's function by 1. Quartic interaction is 
similar, but the process will increase the order by 2. 

(2) Add imaginary unit i times a sum of the n — 1 
graphs formed by pairing each leg with leg 1 and connect- 
ing with the propagator Go, multiplied by a (n — 2) order 
remaining Green's function. 

(3) Symmetrize the graphs, if desired, i.e., do steps (1) 
and (2) for every leg, 1, 2, • • •, n, add them up, and then 
divide by n. 

These rules can be programmed with a symbolic lan- 
guage such as Mathematica. The results are identical 



to the usual expansion of exp[— i J V n (T)dr] and applying 
the Wick theorem. 



3.2.5 Heat current and conductance 

The formulism discussed above offers us computational 
method for the Green's functions, but the most impor- 
tant quantity to calculate in thermal transport is the heat 
current. In this subsection, we discuss the relationship 
between the Green's functions and heat current or ther- 
mal conductance. The derivation of heat current formula 
is similar to that of electron current in [761128] , For in- 
teracting phonon systems, this has been carried out in 
[121131130183] . 

The average rate of energy decrease in the left lead is 
the current flow from the left lead to the central region, 



It. 



(H L (t))- 



(131) 



In steady state, energy conservation means that Jx, + Ir = 
0. Using the Heisenberg equation of motion, we obtain, at 
t = 0, I L = ((u L ) T V LC u c ). The expectation value can be 
expressed in terms of the Green's function G^ L (i,t') = 
—i{u L (t')u c {t) T ) T . Since operators u and u are related in 
Fourier space as u[lu] = — iwu[w], we can eliminate the 
derivative and get, 



1 f°° 

— y Tr(V LC G< L H) W ^. 



(132) 



Using the result derived earlier relating the mixed lead- 
center Green's function to the center-only Green's func- 
tion (transpose of Eq. ([12411 with Go replaced by G) and 
applying the Langreth theorem, Eq. we have G^ L [oj] = 
G r cc [uj}V CL g^[uj} + G^ c [uj}V CL gl[uj}. The expression for 
the energy current is then given by 



II = - 



1 f +oc t 

- dcuuj Tr( G r [uj] Ef [uj]+G < [w] Y, a L [u] 



(133) 



where Zl(t, t') = V cl q l {t, t')V lc . For notational sim- 
plicity, we have dropped the subscript G on the Green's 
functions denoting the central region. We can obtain a 
symmetrized expression with respect to left and right lead 
and make it explicitly real, 



dtu 



{(G r - G a )(S< - S<) + iG<(r R -T L )}, (134) 



where r a = i{S r a - ££). Equation [[T33]) or (fT34]) is valid 
for any finite temperature differences of the leads. 
We define the thermal conductance as 



r 1 61 
a = lim — — = — , 

at^o AT ST 



(135) 



where AT is the difference of the temperatures between 
the leads, such that T L = T + AT/2 and T R = T- AT/2. 
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Fig. 7. The 13th Feynman diagram for the nonlinear self- 
energy S n with labelling. 



To take the limit AT — > explicitly, we introduce vari- 
ational derivatives of the Green's functions, self-energies, 
or current, in the sense, e.g., 



5G 
~5T 



lim 

AT-*0 



G(T L ,T R ) - G(T,T) 
Tl-T r 



(136) 



This is different from varying the temperature of the whole 
system. The variational derivation measures the response 
of the system with a fixed overall temperature T and vary- 
ing only the temperatures of the leads with respective to 
an average temperature. With the variation derivative de- 
fined, we can express the conductance in a form similar 
to the Landauer formula ([S]) for the ballistic transport as 
PIT36], 

i r, ^ 



2tt 



du) loT[lu]- 



dT 



(137) 



with an effective transmission coefficient, 

f[ui\ = ^T r {G r (r L + ir n - s)G a r R } + 

^TY{G a r L G r (r R + ±r n + s)}, (138) 

where the nonlinear effect is reflected in the extra terms, 
r n = i{S r n - X7«) and 



S = 



f sr n 5E< 
1 ST 1 5T 



dT 



(139) 



When T n = S = 0, we recover the Caroli formula. We 
note that the effective transmission T[u>\ is temperature 
dependent, while the ballistic transmission is not. 



3.2.6 Perturbative and mean-field self-energies 

The nonlinear self-energy, £ n (r, r'), is the key quantity in 
solving the heat transport problem in interacting systems. 
Clearly, it is impossible to evaluate it exactly. Some form 
of approximation is unavoidable. As the nonlinear inter- 
actions are relatively weak for real systems, perturbative 
treatment can be reliable for temperatures not too high. 
The approximation can be improved, at least in principle, 
if more diagrams can be included. 

The diagrams shown in Fig. [5] are for the contour- 
ordered functions. The rules of Feynman diagrams to for- 
mula are as follows: label each leg of a vertex with (ji,n) 
where / runs over the legs of all the vertices. A vertex is 
associated with Tj t j j (77, ry , 77") for a three-body in- 
teraction, and similarly for a four-body interaction. If leg 




400 600 

Temperature (K) 

Fig. 8. Thermal conductance as a function of temperature for 
a ID junction with three atoms in different approximations. 
The harmonic spring constants are fct = 1.56, kn — 1.44, 
kc = 1.38 (eV/(A 2 amu)). The coupling between the center 
and leads is the same as the spring constants of the respec- 
tive leads. The nonlinear strength is t = 1.8 eV/(A 3 amu 3 / 2 ), 
with an interaction of the form (t/3) ^(m, — iij+1) 3 . Small on- 
site quadratic potentials are applied to the leads with spring 
constants k°™ itc = 0.006, and fc™ sito = 0.013 (eV/(A 2 amu)). 
From ref. |12| . The squares and crosses are the molecular dy- 
namics results for the ballistic and nonlinear systems, respec- 
tively, from this work. The MD step size is At = 3 x 10~ 16 s 
with 5 x 10 8 MD steps. 

I and I' are connected by a line, associate them with the 
Green's function G° ^ (77, Ty). Integrate on the contour 
all the internal r variables and sum over the internal site 
indices j. Since each n-leg vertex has n — 1 generalized 
delta functions in r, this effectively makes each vertex 
with one single r. For example, the 13-th diagram for the 
self-energy is labelled as shown in Fig. \7\ then the associ- 
ated expression is 

J ^ r l^ji ,33 J4 Js^ie ,h Js^jg Jio ,32^j 3 ,j 10 ( T ' T ) 

G J Bj6 (^n)G° 4j7 (r,r 1 )G« 8j9 (r 1 ,r'). (140) 

High order perturbative treatment is computationally 
too costly except for very small systems. Mean-field ap- 
proximations are usually used. In such approximations, 
only certain selected diagrams are used, with the replace- 
ment of G° by the full Green's function G. So the self- 
energy is considered a functional of G instead G°. For 
Coulomb-type interactions in electronic transport, GW or 
even simpler treatments have been used [1461147] . Certain 
types of self-consistent approximation have the added ad- 
vantage that the energy conservation [1481145] is fulfilled 
in the sense II + Ir = 0. This may not be the case if bare 
self- energies are used. 

In ref. [130j , only the first self-energy diagram in Fig. [5] 
is used. We have made a comparison of the first and sec- 
ond order perturbative approximation, and self-consistent 
mean- field approximation (with the first two diagrams), 



J3 'J4 )J5 'J6 
J7'J8'J9>J10 



Jian-Sheng Wang et al.: Quantum thermal transport in nanostructures 



17 



see Fig. [H Our experience with several models [12 13 
seems to suggest that mean-field approximation tends to 
over estimate the nonlinear effect. Additional problem as- 
sociated with a self-consistent approximation is that the 
iterations may not converge. Thus, finding a good approx- 
imation is still an open challenge for the phonon interac- 
tions at high temperatures. 



3.3 Classical molecular dynamics with quantum 
heat-baths 

Molecular dynamics (MD) has been used extensively to 
study thermal transport based on Green-Kubo formula 
in equilibrium or for noncquilibrium MD with heat source 
and sink (heat baths). However, since MD follows a purely 
classical dynamics, quantum effects are not taken into ac- 
count. There have been some attempts to partially incor- 
porate quantum feature into an MD by scaling the temper- 
ature and the thermal conductivity after the simulations 
|149|150|151j . Given the actual classical simulation tem- 
perature Tmd with kinetic energy ./V/cbTmd/2, where N 
is the number of degrees of freedom, a new temperature 
T can be defined by comparing to a quantum harmonic 
system with the same kinetic energy, 



Nk B T } 



MD 



i 



exp(huJk/(kBT)) - 1 



where the summation is over the vibrational modes. The 
conductivity is also scaled by |150j 



K — KMD - 



MD 



dT 



(142) 



which is effectively multiplying the classical value by the 
quantum heat capacity of the corresponding harmonic lat- 
tice. Both Eq. (|14ip and (|142j) are plausible but lack a 
more rigorous basis. Another approach is the centroid 
molecular dynamics proposed by Cao and Voth [152|153j . 
This appears very sophisticated. The main difficulty here 
is to obtain the temperature-dependent effective potential 
for the molecular dynamics. 

In ref . [154] , we proposed a method based on a general- 
ized Langevin dynamics. The quantum Langevin dynam- 
ics is well studied |155ll56j . However, its use in thermal 
transport is only recent [157 82J . The idea is to consider 
a quantum Langevin equation where the degrees of free- 
dom of the leads are eliminated in favoring that of only 
the central interacting region. This can be performed ex- 
actly if the leads and the couplings between the center and 
leads are linear. Then the central part is treated as a clas- 
sical system, while the leads retain quantum-mechanical 
properties. This amounts to a quasi-classical approxima- 
tion of the quantum Langevin equation [158, 159J . Such a 
treatment also works for electron transport, and electron- 
phonon interactions |160j . 

With the Hamiltonian defined by Eq. ([!]), the corre- 
sponding quasi-classical generalized Langevin equation is 




Fig. 9. Energy (a) and heat capacity (b) per degree of freedom 
for the ID chain models studied in ref. 154 . The solid lines are 
exact results for a ID harmonic chain of length N = 128 with 
periodic boundary condition, while the circles are obtained 
for the same model from the MD with quantum heat-baths 
at the ends; the crosses are the MD results of a nonlinear 

with 



model with an extra quartic onsite potential, 
H = leV/(A 4 amu 2 ). 

given by [154] 



u c = -K c u c + F n (u c )- S r (t,t')u c (t')dt' + Z l +Zr, 

J to 

(143) 

where F n = —dV n /du is the nonlinear force, S r is the 
retarded self-energy of the leads, S r = S r L + E r Rl as used 
in the NEGF calculation, but in the time domain; S r L — 
V CL g r L V LC . A similar equation holds for the right lead 
S R using the right lead surface Green's function. 

The noise £z,(f) or £r(£) has zero mean and a cor- 
relation matrix determined by the temperature and the 
model of the leads. Using the (surface) density of states, 
the expression for the correlation can be simplified to get 
a rather compact result for the spectrum of the noises [82] , 



f L (u) + \)hr L [u>], 

(144) 

where r L [oj] = i(E r L [u) - £%[u]) = -2ImV CL g r L {uj}V LC , 
and /l(w) is the Bose- Einstein distribution function at the 
temperature of left lead. The situation for the right lead is 
analogous. The spectrum function F[u>] is even in w and 
is a symmetric, positive semidefinite matrix. A classical 
limit is obtained if we take + 1/2) 7i » fesT^/a;. 

The steady state thermal current can be computed in 
several equivalent ways: 

II = -In = -<^> = ((u L fV LC u c ) 

= -(u c (t) T B(t)) = {u C (t) T B(t)), (145) 

where B(t) = - j£ 27£(t, t')u c '(t')dt 1 + &(t). 

If V n = the linear Langevin equation can be solved 
analytically. The result for the thermal current reproduces 
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exactly the Caroli formula [82] • On the other hand, if tem- 
perature is sufficiently high (e.g., 1/2 of the Debye temper- 
ature), the nonlinear Langcvin equation agrees with clas- 
sical MD result (with appropriate heat baths). For non- 
linear systems with V n ^ 0, the equation can be solved 
numerically. At low temperatures, the transport is mostly 
ballistic, the dynamics correctly reproduces quantum re- 
sults and agrees well with NEGF calculation [154 . At 
the intermediate temperature range, where both quan- 
tum and nonlinear effects are strong, the dynamics gives 
only approximate results. However, we note that the quasi- 
classical approximation can be systematically improved, if 



the form, dP^/dt = 7(P /J _i — 2P [i + P^+i), can be de- 
rived, where P M is the probability that site /i is at ex- 
cited states. The diffusion constant 7 is then related to 
the thermal conductivity. Segal et al. studied thermal rec- 
tifying effects of quantum systems with somewhat similar 
approach of subspace projection, where again diffusion- 
like rate equations are derived for the occupation of the 
quantum states [ 1 7311 741 1 75| . the nonlinearity is intro- 
duced by restricting to two-level or finite-level systems. 
In ref. [176| a master equation approach for more gen- 
eral systems is considered and the Fourier's law of heat 
conduction is derived. The approximations made always 



we use the techniques proposed in [T61|162|163|164|165|166| lfel tOftjc systems in diffusive regime. It seems to us that 



For example, the leading order correction is obtained if we 
enlarge the equation set to include variables u 2 , up, and 
p 2 (p — ii) as independent variables. Thus, the dynamics 
can provide a nonperturbative solution to the quantum 
transport problems. 

In Fig. [8] we show comparisons of MD and NEGF re- 
sults. For the ballistic case, the agreement is perfect with 
very high precision. For the nonlinear system, in order to 
have a stable MD, we have added an additional quartic 
potential, (J-/4:^2(uj ~~ u j+i) i -> wltn A 4 = 0.6 eV/(A 4 amu 2 ) 
in the simulation, thus the comparison is hardly on an 
equal footing. The discrepancy at low temperatures be- 
tween NEGF and the MD has to do with the zero-point 
motion contribution. In the noise spectrum, Eq. (|144[) . 
the 1/2 in Jl + 1/2 represents a contribution from zero- 
point vibration of the quantum oscillators. We can replace 
/l + 1/2 by /l(w) for ui > and — u) for uj < 0, and 
it will not have any effect to the thermal current of a bal- 
listic system. But the removal of the zero-point motions 
does have some consequences to the nonlinear system un- 
der quasi-class approximation. 

It is interesting to note that the Langevin dynamics 
with quantum heat baths also provides a way to compute 
the equilibrium energy, and by numerical differentiation, 
the heat capacity of the quantum models. The energy is 
obtained by the MD averages of the classical expression 
of kinetic energy plus the potential energy, with a left 
and right bath at the same temperature. A comparison 
of the exact and MD results is presented in Fig. [9] for the 
linear chain. For nonlinear systems in thermal equilibrium, 
an alternative method is (path-integral) quantum Monte 
Carlo [T69] , 



3.4 Other treatments of nonlinear effects 

In this subsection, we briefly mention several other ap- 
proaches for treating the nonlinear quantum thermal trans- 
port for dielectric materials. In refs. [17011711172) Michel 
et al. used Hilbert space average method to derive diffusion- 
like equations on ID chains, thus justified the Fourier's 
law and computed the thermal conductivity from Hamil- 
tonian model parameters. These models are idealized so 
that each site has a band of excited states separated with 
an energy gap to the ground state. The couplings be- 
tween the sites are assumed weak. Diffusion equation of 



it is difficult to address the issue of transition from bal- 
listic behavior at short time and length scales to diffu- 
sive behavior at long time and length scales in the above 
approaches, but see ref. |177j for such possibilities. In 
refs. [178|179|180j heat flow through model molecular and 
nanocrystal junctions is treated using Fermion-goldcn rule 
for the rate of phonon relaxations. The Fermi-Pasta-Ulam 
chain is treated quantum-mechanically in ref. [181] . The 
authors find an L 0A divergence for the thermal conduc- 
tivity by considering the relaxation rates of the phonon 
modes. 



4 Electron-phonon interactions 

In this section, we extend the NEGF method in subsec.[ 
to include the electron subsystem and electron-phonon 
interactions (EPIs). The study of EPI in the context of 
nanoscale electronic transport is an important field which 
has attracted intense research recently [61182] . An exten- 
sive discussion of this issue is certainly out of the scope of 
current review. Here we only concentrate on one special 
aspect, namely the current induced heating in nanostruc- 
tures. Different schemes exist in the literature to study this 
effect |183I184I185I186I187I188I189I19()I191I192I193I194] . When 
the EPI is strong, a canonical transformation is useful to 
study a minimum model system [133 6 83 193J. But it is 
not applicable to large systems. On the other hand, per- 
turbative approaches based on the lowest order Born ap- 
proximation are applicable for relatively weak interactions 
|184ll92l83ll94j . Actually, it has already been used in the 
first principles study [19011911192] . Hybrid approaches also 
exist, where the electrons are treated quantum-mechanically, 
while the phonons use MD with quantum corrections [183118711881189] . 
In our NEGF formalism, we use the mean field method 
based on the lowest order Born approximation, the so- 
called self-consistent Born approximation (SCBA). 

The Hamiltonian of the whole system is the sum of 
the electron and phonon part. Consisting with Eq. |T]), 
the electron Hamiltonian reads 



CC 



a=L,C,R 



a=LM 



ijk 



k u k c]c 



(146) 

where c a is a column vector consisting of all the annihi- 
lation operators in a region, <v a is the corresponding cre- 
ation operators. V" c has similar meaning as V aC in the 
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phonon Hamiltonian, and V" c = (V^ a )^ . The last term 
represents the EPI, where we only keep the first order Tay- 
lor expansion about the atomic equilibrium position. 
is the interaction coefficient. We assume that EPI takes 
place only at the center region, and the subscript C has 
been omitted. 

In ref. [194] we adopt the SCBA to study the coupled 
electron and phonon transport. All the electron Green's 
functions are defined as usual |128J, e.g., 



D r (t,t') = -id(t-t')([c(t),cHt')] + ). 



(147) 



We use d to denote the free lead Green's function. Do and 
D are the noninteracting and the full electron Green's 
functions of the center, respectively. The free lead Green's 
function d can be obtained using the same method as that 
for phonons (a = L, R) 



d<[w] = -fl[u](f a 



uA-dtU 



(148) 
(149) 



f e (j-o) is the Fcrmi-Dirac distribution function. The non- 
interacting electron Green's functions can also be solved 
exactly 

Dl' a [w]=[(u±i V )-F°-ir> a [u,]]- 1 , (150) 

D<<>[u,]=ir [w]ii < ' > [u]r%H (151) 

The electron self-energy includes contributions from both 
leads II — IIl + U n, and 

n^<>[u>] = v? a d r f'<>[u]v: c . 



(152) 



The full Green's function D follows the same relations, 
except that the self-energies include EPI contribution. The 
notations used in Eq. (|147[) to (| 1 52[) for electron Green's 
functions and self-energies are not standard. As we have 
already used G and £ for the phonons in earlier sections, 
D and 77 are used here for the electrons. 

The central problem is how to get 77 op i and S ep i- The 
SCBA includes only the lowest order self-energies. For 
iTcpi, there are two Feynman diagrams, the first two in 
Fig. [5] with modifications. In the first one, only half of 
the circle represents phonons, and all other lines repre- 
sent electrons. In the second one, only the vertical line 
represents phonons, and all others are electrons. For i7 ep i, 
there is only one diagram, the first one in Fig. [5] with the 
circle representing electrons. The coefficients also change 
to i, —i, and — i, respectively. The rules to translate the 
diagrams into expressions remain the same as for phonons, 
provided that the interaction vertex is defined by 



M t k Jn,T 3 ,T k ) = Mh 5 (t„- , Tj ) <5 (t,- , Tfe ) 



(153) 



Mean field method is used to calculate the self-energies 
and the Green's functions. 

The electron and the energy current can be expressed 
via the Green's functions. The electronic current out of 
the lead a (a = L, R) is [7611271128] 

/+oo J 
^ —Ti{D > [uj}n<[u}}-D < [uj)n>[uj]}. (154) 



Similarly, the energy current is 
" + °° dcu 



2vr 



Tr{D> M7T< [u] - D< [u)II> [w] } . (155) 



The heat current can be obtained from Eqs. (|154H155| as 
la = 1% ~ HaJa/e- Ha is the chemical potential of the lead 
a. For phonons the energy and the heat current are the 
same. In the presence of EPI, there is energy exchange 
between the electron and the phonon subsystem, which is 
the source of Joule heating. The use of SCBA guarantees 
the energy and electronic current conservation [1941147) . 
so that we can write them in symmetric forms. The elec- 
tronic current is 



+ oo 



^T>](/£M-/^M). (156) 

The effective transmission coefficient reads (u) omitted) 
f e = Tr{±[D r (Al + ±A cpi - X cpi )D a A c R 

+D a A c L D r (A c R + ±A cpi + X cpi )]}, (157) 

5(/* + /EMc P i + ^ c < pi 



(158) 



fe fe 

J L JR 

Al = z(n r a - n a a ) (a = L,R) and yl opi = z(77^ pi - 77 c a pi ). 
The total energy current includes electron and phonon 
contributions 



|^{f>](/£M-/£H) 
+ifp h H(/ i H-/ i? H)}, (159) 



where T ph is similar to eq. (|138[) . except that F n is re- 
placed by Tcpi = i{S r cpi - Z'epi)' and S b y 



UfR + fL)r epi -iZ< pi 

II - fn 



(160) 



Here we have corrected a sign error in Eq. (22) of ref. [194] . 
Heat generation during the current flow can also be writ- 
ten in terms of Green's functions [194] . 



jn- 



-oo -oo nmikl 

(161) 

Equation (1161|) can be used to study the Joule heating in 
noncquilibrium molecular junctions. Similar results have 
also been derived by other authors 83 192 193 . 

In Fig. [TU1 we show the heat generation as a function 
of onsite energy of a single quantum dot coupled with two 
ID leads, calculated using Eq. (|161|) . When the applied 
bias is less than the phonon energy to, there is no heating. 
When the bias energy is slightly larger than u> and less 
than 2lu, there are two energies where the heat generation 
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e(eV) 

Fig. 10. Heat generation Q as a function of electrical onsite 
energy e for a quantum dot coupled with ID leads (ref. [194] ). 
From the inner to the outer side, the applied biases are V — 
0.05 to 0.50 V with an increment of 0.05 V, respectively. 

is the largest. These two peak positions are approximately 
at —0.5 eV + u> and 0.5 eV — u>. They merge into a single 
one at a bias of 2u> until it reaches saturation. After that, 
this peak broadens to a ladder. These interesting features 
are unique properties of nanostructures, absent in the bulk 
case. 



5 Experiments in nanostructure heat 
transport 

Some of the earliest direct mesoscopic thermal conduc- 
tance measurements are done on patterned semiconductor 
structures with cross section sizes of few hundred nanome- 
ters 195 . The calorimetry on such a small scale is a chal- 
lenge, but very promising, e.g., one may count the number 
of phonons in a nanostructure 196 1. In 2000, such type of 
devices is used to confirm the universal thermal conduc- 
tance quanta at temperatures below 1 K . Bridges of 
catenoidal shape are used, based on an earlier theoretical 
analysis that such geometric shape gives the maximum 
transmission [17] . Other extensively studied systems are 
carbon nanotubes. Initially, only mats or rope bundle or 
film form can be measured, given a relatively low ther- 
mal conductivity of order 100 W/(mK) at room tempera- 
ture [197 198 199J. The initial theoretical predictions by 
molecular dynamics for single nanotubes along the ax- 
ial direction are much higher 200 201 1. Indeed, individ- 
ual multiwalled and single wall carbon nanotube measure- 
ments became possible, and the measured values are about 
3000 W/(mK) at room temperature by several groups 
[202 120312041205] . Typically, these measurements use tube 
lengths of few j«m. The length dependence effects are an- 
alyzed in [205] and observed in [206] . In ref. [207] . it is 
claimed that ballistic phonon transport is consistent with 
experimental data, even for temperatures as high as 900 K 
for submicron tube lengths. Another class of systems is 



individual silicon nanowires 208 209 210 . It has been ob- 
served that rough and small diameter silicon nanowires 
have a low thermal conductivity; this property can be used 
for good thermoelectric applications [2T0] . Thermal recti- 
fication effect has been observed in carbon nanotubes with 
one side of the tube deposited with CgHi 6 Pt [211] , How- 
ever, the observed effect is rather small, at a level of few 
percents. Another way of explicit control of the heat flow 
in nanostructure is using a tunable thermal links 212J, 
where a multiwalled carbon nanotube is cut and the outer 
layers of the tube can be moved mechanically with re- 
spect to the inner layers. Much large modulation is then 
possible. In ref. [213] . it is demonstrated that the carbon 
nanotube heat conduction is rather robust again bending 
and deformation. Phonon transport in point contact of 
metal-insulator is studied in ref. [214] . The studies of ther- 
mal transport at an individual molecular level are scarce. 
But impressive progress is made recently [21512161217] , 
In ref. [217] . alkane chain anchored to a gold substrate is 
heated to about 1100 K, the propagation of the heat waves 
is timed. This timing information indicates that the heat 
transport at such scale of few nanometers is ballistic. The 
thermal conductance was inferred to be 0.05 nW/K. 

6 Conclusion 

Using a general junction model as an illustration, we con- 
sider the problem of theoretical modeling of heat transport 
in nanojunctions. If the nonlinear forces can be omitted, 
this would be the case if temperature is low or the sizes 
of the system is small, the Landauer formula gives a very 
good description of the ballistic thermal transport. We 
give a derivation of the Landauer formula from the wave 
scattering point of view. A derivation from NEGF is also 
outlined as a special case of the more general nonlinear 
result. Several different ways of calculating the transmis- 
sion coefficient appearing in the Landauer formula are re- 
viewed. Among the alternative computational methods, 
the calculation based on Caroli formula and the iterative 
algorithm of the surface Green's functions is the most ef- 
ficient computationally for the total transmission coeffi- 
cient T[u>]. The mode-matching method gives the results 
of transmission of each mode to every other mode, thus 
containing more information. The relation of the two ap- 
proaches is clarified. The transfer matrix method and Pic- 
card formula are theoretically interesting, but they may 
not be practical in actual calculations, e.g., the transfer 
matrix can be numerically unstable for large systems. In 
summary, for the ballistic transport, efficient computa- 
tional algorithms exist for arbitrary complex systems with 
numbers of atoms less than about 10 3 . 

For the treatment of nonlinear effects, it is still a great 
challenge. We first presented a phenomenological treat- 
ment of the nonlinear effect using the information on the 
mean free path. A systematic development of the nonequi- 
librium Green's function method is given. We emphasize 
the use of contour ordered Green's functions and the equa- 
tions of motion defined on the Keldysh contour. To do 
this, we must generalize the 9 and S functions, and also 
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generalize the Heisenberg equation to the contour. Quan- 
tum field-theoretic method is used to develop the Dyson 
equations relating the ballistic system and interacting sys- 
tem. The nonlinear effect can be taken into account if 
the nonlinear self energy S n can be obtained accurately. 
Perturbative expansion and mean-field approximation to 
the nonlinear self-energy are possible. The Landauer for- 
mula can be generalized formally in the nonlinear case, 
except now the transmission function is temperature de- 
pendent. For small systems at not too high temperatures 
(e.g., T < 400 K), the results with these approximations 
are reliable. But for large systems or higher temperatures, 
the approximation can break down. A more promising 
method to handle large systems and high temperatures 
is molecular dynamics. With the quantum heat baths dis- 
cussed in this review, the quantum effect can be partially 
taken into account and it is even possible to correct the 
quasi-classical approximation. 

We also reviewed briefly methods for electron-phonon 
interactions. This is an important area where many inter- 
esting effects can be investigated experimentally, and it is 
also very relevant to the electronic industry. NEGF and 
other methods will play important role in the study. We 
also give a quick review of the experimental status of ther- 
mal transport problems in nanostructures. To a theoreti- 
cal or computational modelling scientist, it is easy to work 
with small systems, while it is hard to an experimentalist 
to manipulate and measure heat at the molecular level. 
We hope that the two sides of approaches will converge 
well in the near future. 

For the convenience of reader, we have put up a web- 
page |218j to post some of the implementations of algo- 
rithms, which may clarify the doubts and subtlties in un- 
derstanding them. 
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